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Abstract 


The  eigenvalue  method,  as  presented  by  T.  M.  Shlh  and 
J.  T.  Skladany,  is  evaluated  to  determine  the  advantages  and 
disadvantages  of  the  method  as  compared  to  fully  explicit, 
fully  implicit,  and  Crank-Nlcolson  methods.  Time 
comparisons  and  accuracy  comparisons  are  made  in  an  effort 
to  rank  the  eigenvalue  method  in  relation  to  the  comparison 
schemes . 

Shlh  and  Skladany's  original  results  are  verified  by 
duplicating  their  efforts  with  the  method.  The  eigenvalue 
method  is  used  to  solve  the  parablolc  heat  equation  in 
multidimensions  with  transient  temperatures.  Extensions 
into  three  dimensions  are  made  to  determine  the  method's 
feasibility  in  handling  large  geometry  problems  requiring 
great  numbers  of  internal  mesh  points. 

The  eigenvalue  method  proves  to  be  slightly  better  in 
accuracy  than  the  comparison  routines  because  of  an  exact 
treatment,  as  opposed  to  a  numerical  approximation,  of  the 
time  derivative  in  the  heat  equation.  It  is  an 
unconditionally  stable  method.  It  has  the  potential  of 
being  a  very  powerful  routine  in  solving  long  transient  type 
problems.  The  method  is  not  well  suited  to  finely  meshed 
grid  arrays  or  large  regions  because  of  the  time  and  memory 
requirements  necessary  for  calculating  large  sets  of 
eigenvalues  and  eigenvectors. 


EVALUATION  OF  THE  EIGENVALUE  METHOD  IN 

THE  SOLUTION  OF  TRANSIENT  HEAT  CONDUCTION  PROBLEMS 

I.  INTRODUCTION 

Heat  flow  in  solid  bodies  and  across  their  boundaries 
is  a  common  engineering  problem.  The  solution  of  these 
problems  can  be  found  by  direct  calculation,  experiment,  or 
by  use  of  numerical  techniques.  Calculation  is  the  most 
preferred  method,  but  analytical  solutions  may  not  exist  in 
many  problems.  Experimentation  can  be  expensive  and  is 
often  inconvenient,  if  not  impossible.  Numerical  analysis 
becomes  a  very  good  alternative  to  the  problem  solver. 

At  first  glance,  numerical  analysis  of  heat  flow  seems 
"too  complicated"  and  "can't  be  trusted,  as  it  is  only  an 
approximation  anyway"  (9:viii).  This  is  not  a  fair  opinion 
as  there  are  many  useful,  simple,  and  very  accurate 
techniques  available.  Among  these  are  the  fully  explicit 
(7:223-251),  fully  implicit  (1:49-53),  Crank -Nicol son 
(8:50-67),  and,  a  fairly  new  thechnique  in  heat  transfer, 
the  eigenvalue  method  (19:409).  These  do  not  comprise  a 
complete  list  of  techniques  that  are  available,  but,  rather 
represent  broad  classes  of  approaches  in  discretized 
solutions. 

The  eigenvalue  method,  "which  has  been  used  by 
researchers  in  structure  mechanics,  but  is  relatively 
unpopular  in  the  heat  transfer  community"  (19:409) ,  can. 


indeed/  be  applied  to  heat  transfer.  It  is  especially  well 
suited  to  heat  transfer  problems  with  long  transient  times 
until  equilibrium  is  reached.  It  differs  from  the  other 
techniques  in  that  eigenvalues  and  eigenvectors  must  be 
calculated.  Once  this  hurdle#  albeit  a  substantial  one#  is 
cleared#  the  eigenvalue  method  uses  algebraic  type  equations 
to  arrive  at  its  solutions#  thus  avoiding  time  consuming 
array  solving  or  iterations  to  converge  to  correct  answers. 

PROBLEM 

A  recent  (1983)  article  published  by  Shih  and  Skladany 
(19:409-422)  investigated  the  eigenvalue  method  to  solve 
transient  heat  flow  problems.  The  problem  Investigated  in 
this  study  is  Independent  verification  of  their  results. 
Additionally#  extensions  of  the  eigenvalue  method  will  be 
made  into  three  dimensions  and  the  method  will  be  used  in  an 
Independent  problem.  Advantages  and  disadvantages  of  the 
method  will  be  identified.  The  method  will  be  compared  to 
the  fully  implicit#  fully  explicit#  and  Crank-Nicolson 
methods  and#  finally,  an  overall  rating  will  be  given  In 
relation  to  these  methods. 

SCOPE 

This  study  will  use  cartesian  coordinates  in  two  and 
three  dimensions.  In  addition  to  the  comparison  routines# 
each  method's  solution  will  be  compared  to  the  exact 
analytical  solution  using  a  least  squares  analysis.  The 


eigenvalue  method  will  be  extended  to  three  dimensions  and 
compared  to  a  three  dimensions  version  of  the  fully  implicit 
method. 

Cylindrical  or  spherical  coordinates  will  not  be 
Incorporated  In  the  study.  It  is  not  the  purpose  of  this 
study  to  perform  a  detailed  analysis  of  the  comparison 
routines.  All  problems  investigated  will  consider  only 
Dlrlchlet  boundary  conditions  to  Insure  uniform  comparisons. 
Problems  Involving  heat  sources  within  the  solid  will  not  be 
examined . 

GENERAL  APPROACH 

The  first  objective  will  be  to  explain/  in  detail/  what 
the  eigenvalue  method  (EM)  Is  and  why  It  works.  Next/  the 
exact  analytical  solution  to  the  two  dimensional  problem 
presented  In  the  Shih  and  Skladany  article  will  be  verified 
(19:410) . 

Armed  with  this  Information/  the  results  of  the  article 
will  be  recomputed  based  upon  the  sample  problem  they 
prepared.  The  same  problem  will  be  solved  using  the  three 
comparison  techniques:  fully  explicit  (EX)/  fully  implicit 
(IM) /  and  Crank-Nlcolson  (CN) . 

Each  run  will  be  timed  for  wall  clock  times /  user  time 
(time  actually  obeying  commands) /  and  central  processor  unit 
time  (CPU  time)  (16:473-474).  Each  method's  output  will  be 
compared  against  the  exact  solution  and  a  least  squares 
analysis  will  be  performed  to  produce  an  overall  accuracy 


for  each  run  (23tl021-'1035) .  All  runs  will  be  made  using 
the  AFIT  VAX  11/780  computer  system. 

The  sample  problem  will  be  extended  into  three 
dimensions  to  investigate  ease-o£-use  of  the  EM  in 
multidimensions  and  perform  an  problem  independent  from  that 
of  the  article.  An  additional  three  dimensional  problem 
will  be  performed  to  deliberately  remove  any  symmetry  for 
the  purpose  of  evalutatlng  the  EM  in  these  type  cases. 

Finally,  the  EM  will  be  closely  examined  for 
enhancements,  strengths,  weaknesses,  etc.  An  overall 
performance  rating,  in  relation  to  the  comparison  routines, 
will  be  given. 


II.  EXISTING  THEORY 


Heat  conduction  in  solids  has  been  the  object  of 
investigation  by  engineers  for  quite  some  time.  The  heat 
diffusion  equation  is  given  (21:9-54)  as 


+  (1) 


in  three  dimensions  where 

^  >  material  density 

Cv  *  specific  heat 
^  «  thermal  condictivity 

All  problems  examined  in  this  study  will  be  based  upon 
this  equation.  Shih  and  Skladany  used  a  two  dimensional 
form  of  Eq.  (1)  in  their  investigations  (19:410-413).  For 
simplicity,  the  following  three  comparison  routines  will  be 
examined  in  just  the  two  dimensions. 

EXPLICIT 

One  of  the  simpler  numerical  techniques  is  the  explicit 
formulation.  The  relative  ease  of  setting  the  explicit 
method  up  and  the  absence  of  large  simultaneous  equations 
are  clear  advantages  of  this  method  (5:91).  Its  major 
disadvantage  lies  in  the  commonly  know  problem  of  stability 
(19:409;  7:377).  This  problem  links  the  time  step  of  the 


■ethod  to  the  spatial  step.  This  can  prove  to  be 
prohibitive  in  large  transient  type  problems. 

Use  of  standard  finite  differencing  techniques  will  be 
used  in  all  formulations  in  this  study  (7:376-378;  5:59-63) 
Assuming  a  forward  difference  approximation  for  the  time 
derivative  and  a  central  difference  approximation  for  the 
spatial  derivatives,  the  explicit  method  can  be  written  as 

/Ajf  +  O(hi^)  (4) 

and  plugging  into  the  heat  equation  (1) 

^TA-l  »  (X  ,5, 

where 

=  b/  fCv 


yields  the  explicit  method  algorithm.  Placing  all  the 
knowns  on  the  right  side,  the  remaining  unknown  is  the 
temperature  at  each  node  for  the  new  time  step. 


Let^  ■Ay  ■  h,  then 


+  t4.  -  ^Tijl  +  Tifj  (6) 

This  algorithm  can  be  used  to  calculate  temperatures  at 
any  nodal  point  (l,j)  at  any  time  t,  provided  the  At  Is 
sufficiently  small  to  meet  the  stability  requirement. 

To  Insure  a  bounded  answer »  but  not  necessarily 
accuracy  (15),  the  condition 

At/  ^  I  A| 

must  be  met  (7:379**380) .  This  method  has  truncation  errors 
with  order  0(ht.)  >  0(h|^)  (5:59-63). 

IMPLICIT 

To  avoid  the  serious  stability  criteria  of  the  explicit 
method,  the  engineer  can  turn  to  an  implicit  formulation. 

The  main  principle  of  this  scheme  is  that  the  temperature  at 
each  point  Is  found  at  all  nodes  simultaneously  using  the 
temperature  values  from  the  previous  time  step  (17:58). 

This  also  gives  rise  to  Its  disadvantage,  in  that  large 
systems  of  simultaneous  equations  must  be  solved  to  arrive 
at  a  solution.  For  finely  meshed  problems  with  long 
transient  times  this  can  become  a  significant  problem. 


Formulation  of  the  Implicit  scheme  is  similar  to  that 
of  the  explicit  scheme  except  the  spatial  derivatives  are 
evaluated  at  the  new  time  step  (evaluated  at  the  current 
time  step  in  explicit) . 

The  time  derivative  will  remain  unchanged  and  the 
spatial  derivatives  will  become  (5:104;  1:49-53) 

CT^j-2Ty'+T^j)/Ax^  +  0((.?)  (8) 

(TiJ'  -2Ty  +t5.)/ai,'  +  Oa?)  (9) 


Plugging  Eqs.  (2)/  (8)»  and  (9)  into  Eq.  (5)  and  then 
rearranging  all  the  unkno%ms  on  one  side  yields 


•4fl 


i*' 


-I 


(10) 


This  method  is  well  adapted  to  the  computer  because  it 
is  unconditionally  stable r  but  it  does  generate  the  large 
sets  of  simultaneous  equations  (11:455).  This  method  has 
truncation  error  on  order  of  0(h4)  +  0(h|^)  (5:59-63). 
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CRANK-NICOLSON 

An  attempt  to  avoid  numerical  instability  led  to  the 
implicit  method.  The  method  is  unconditionally  stable ,  but 
the  method  still  had  truncation  errors  of  0(ht)  O(hM^) 
(5:59-63).  This  can  be  improved  to  0(ht)  *  0(h^ )  by 
approximating  the  time  derivative  with  a  central  difference 
(11:455).  This  method,  it  turns  out,  has  problems  with 
stability  (11:455) .  Crank  and  Nicolson  suggested  to 
overcome  this  by  representing  the  derivative  with  an  average 
of  a  forward  difference  and  a  backward  difference,  each 
weighted  by  1/2  (8:50;  11:455-457).  If 

^  (11) 

then  the  Crank-Nicolson  scheme  is 

T^i'  =  1 +  {f/ C2(i+arn 

This  method  also  generates  large  systems  of 
simultaneous  equations.  The  advantage  comes  from  the 

a 

increase  in  the  accuracy  gained  by  the  0(ht)  term. 


III.  NEW  THEORY 


THE  EIGENVALOE  METHOD 

The  eigenvalue  method  (EM)  is  not  really  a  "new"  theory 
or  discovery.  Shih  and  Skladany  suggest  in  their  article 
that  this  method  has  been  used  extensively  in  structure 
mechanics  (19:405).  The  use  of  eigenvalues  to  solve  systems 
of  equations  has  been  studied  in  math  text  books  for  quite 
some  time.  It  is  the  application  of  the  EM  to  solutions  of 
the  heat  equation  that  is  fairly  new. 


WHAT  IS  T^  ^ 

Solving  problems  numerically  requires  the  division  of 
the  region  in  question  into  small  grids.  The  intersection 
of  the  grid  lines  (called  nodal  points) r  are^  typically,  the 
points  of  Interest  in  the  study.  If  these  grids  are 
sufficiently  small,  the  exact  behavior  under  study  can  be 
quite  accurately  approximated. 

Each  node  point  is  an  unknown  value  and  is  generally 
represented  by  a  governing  equation.  A  system  of  equations 
such  as  these  can  be  solved  simultaneously  to  find  these 
unknowns.  It  takes  n  equations  to  find  n  unknowns. 

Provided  the  equations  can  be  found,  all  that  remains  is  the 
selection  of  a  techlnque  to  solve  them.  Several  techniques 
exist  to  do  just  that;  among  them  is  the  EM. 
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WHY  DOES  THE  EM  WORK? 

For  clarity,  the  EM  is  applied  to  the  transient  heat 
equation.  This  application  demonstrates  the  method  and  why 
it  is  a  valid  numerical  technique. 

The  transient  paraJaolic  heat  equation  in  three 
dimensions  appears  as 

4  4 

Using  standard  central  differences  for  the  spatial 
derivatives  (7:376-378) 

4  (Tc,i,fcti  -  2Ti,},v4Tc,],ti  /Az^)3 

and,  again,  let  Ax  «  ^  «  Az  *  h. 


an  =6 


Figure  1.  Notation  Convention  for  Discretizations. 

Figure  1  depicts  the  notation  convention  used  in  this 
report  for  all  discretization  schemes.  Using  this  notation, 
Eq.  (14)  can  be  written  as 

?Cv,5Tp Ai  -  ( k/\?)  CTe  +  Tw  +  Tn  +  Tf 

+  Is  +  Tb  "  ^  (15) 


where 


p  ^  1,  2,  3f  ...  fn. 


DERIVATION  OP  METHOD 

Borrowing  mainly  from  the  Shih  and  Skladany  report 


(19:409-422),  the  following  is  the  derivation  of  the 
eigenvalue  method.  Equation  (15)  can  be  expressed  in  matrix 


form  as 


[c][t5  .  [KlfTi  +  \f\ 

where  "CC]  is  an  nxn  identity  matrix  multiplied  by  fhCy,  CKl 
is  an  nxn  sparse  conductance  matrix  whose  diagonal  elements 
are  -4k/h"  (-6k/h  in  three  dimensions)  (19:411).  The  vector 
is  of  length  n  and  incorporates  the  boundary  conditions. 

In  general/  all  matrices  have  n  characteristic 
eigenvalues  and  associated  eigenvectors  (14:345).  Most 
discretization  schemes  do  not  make  use  of  these.  Along  with 
the  appropriate  multiplying  constants  determined  from 
boundary  conditions/  these  eigenvalues  and  eigenvectors  can 
be  used  to  solve  the  system  of  simultaneous  equations  from 
which  they  came  (14:824).  This  is  the  basis  of  the  EM  and 
why  it  works. 

Ose  of  some  simple  matrix  algebra  reveals  the  core  of 
the  method.  If  Eg.  (16)  is  pre-mul tipi led  by  ^^3”^  to  yield 
(again/  following  Shih  and  Skladany)  (19:411) 


(17) 


[cTHfr!-  4  [cr{4j 

and  define 


(18) 


to  get 


4  [g] 


Now,  let 


{t\  =  {ml  4  (t,1 


where  is  some  arbitrary  fxinction  of  time  and  is 

independent  of  time  and  represents  the  steady-state 
solution.  The  derivative  of  Eq.  (21)  with  respect  to  time 


gives 


ji]  = 


and  plugging  Eqs.  (21)  and  (22)  into  Eq.  (20) 


and  if 


then  Eq.  (23)  becomes  just 


[ml  =  [i>]W 


(25 


which  is  a  set  of  "first-order/  linear/  ordinary 
differential  equations"  (19:411).  Appendix  B  details  the 
solution  of  Eq.  (25)  as  applied  to  a  system  of  two  equations 
(without  loss  of  generality) .  Based  upon  this  and  assuming 
the  solution  of  Eq.  (25)  to  be 

(26) 

i 

where  |c,*]  is  the  set  of  multiplying  constants  yet  to  be 
determined  and  fAij  are  the  eigenvalues  of  the  system.  It 
turns  out  the  multiplying  constants  are  actually  a  product 
of  the  eigenvectors  and  another  constant  determined  by  the 
boundary  conditions  of  the  particular  problem.  Once  the 
eigenvalues/  eigenvectors/  and  the  multiplying  constants  are 
determined/  the  Eq.  (25)  can  be  solved.  These  results  are 
used  in  Eq.  (21)  to  find  the  temperature  at  each  nodal 
point.  In  matrix  form  the  EM  can  be  written 


[T;  i  ^  [  Ir-i  +  E  ) 


(27) 


where 

i  ®  3  “  1/  2/  3/  ...  /n. 

and  a^  is  the  i-th  multiplying  coefficient  determined  by 
boundary  conditions  and  v^'*  is  the  j-th  component  of  the  i-th 
eigenvector  (see  List  of  Symbols) . 


From  the  previous  derivation  and  the  Shih  and  Skladany 
article  several  advantages  and  disadvantages  seem  possible. 
Among  the  advantages  are 

1.  Better  accuracy;  no  time  truncation  error. 

2.  Efficient  use  with  long  transients. 

3.  Unconditionally  stable  method. 

4.  Rapid  solution;  few  eigenvectors  needed. 

Two  possible  disadvantages  could  be 

1.  Difficulty  in  finding  and  using  good 
eigenvalue /vector  routines. 

2.  Time  consuming  in  fine  mesh  space  problems. 

The  results  of  this  study  will  verify  or  contradict 
these  and  bring  to  light  any  others  not  mentioned. 


IV.  APPLICATION  TO  SPECIFIC  PROBLEMS 

One  of  the  main  purposes  of  this  study  is  to  verify 
results  achieved  by  Shih  and  Skladany  in  their  article. 

They  presented  a  two  dimensional  sample  problem  (19:412). 
For  the  purpose  of  brevity ,  they  only  examined  this  problem 
by  taking  full  advantage  of  the  symmmetry  Involved.  This 
study  does  the  same  and  will  also  examine  the  problem  with 
the  assumption  that  not  all  problems  will  be  symmetric. 

THE  TEST  PROBLEM 

Consider  a  square  slab  of  length f  L,  on  a  side  with 
constant  temperature,  T^)  ,  on  each  of  the  sides.  Further, 
let  the  region  be  governed  by 

•+*  ^28) 

subject  to  the  boundary  conditions 

T(0,y,t)  -  Tp 
T(L,y,t)  =  To 
T(x,0,t)  *  To 
T(x,L,t)  »  Tq 
and  the  Initial  condition 


T(x,y,0)  »  T, 


Pig.  2.  Domain  £or  2-D  transient  heat  conduction  problem 


Figure  2  shows  graphically  the  sample  problem.  Shih 
and  Skladany  maintain  that  the  temperature  at  the  node  point 
in  the  center  is  equal  to  the  temperature  at  Tto  (19:412). 

It  is  not  necessary  to  make  this  approximation  to  assure  h* 
accuracy.  A  more  accurate  formulation  is  to  make  the 
temperature  at  the  center  seperate  and  add  an  additional 
equation  to  the  resulting  system. 

The  analytical  solution  to  Eq.  (28)  can  be  found  in 
Appendix  A  by  removing  the  z  dependence.  The  normalized 
solution  will  then  be 


(29) 


where 


?  (^/nTr)eKpC-n2n^  /t^^SinCnrrKyO 

A»l 


(29a) 


and 


eo 

Bt  (y j't'i  -  2  (H  /m)  eK?  (-tf t-t A?  )S(  n  [  niry  /u^  (29b) 


with 


n  *  i/  3/  5f  ••• 


for  both  01  and  0^  . 

For  simplicity  let  (19:412) 


(30) 


The  reader  is  advised  to  refer  to  the  list  of  symbols 
found  at  the  beginning  of  this  report  in  order  to  keep  track 
of  the  normalized  and  real  values.  For  generalityr  this 
problem  is  done  with  normalized  values. 

In  the  Interest  of  brevityr  only  a  few  of  the  20 
differential  equations  generated  by  Eq  (15)  (2-D  version) 
are  listed 
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2'rz  -  HTT  -f  aTa 


(31) 


^  i  n  / 

“  "Tim  +Ti(»  + 

1 1*1  -  MTii 

(32) 

At 

~  ”Tis  -f  STn  + 

"ITo  -  ^Ti« 

(33) 

=  ETn  + 

-  ^Tiq 

(34) 

®  Tii  +  2Th 

0 

1 

(35) 

These  20  equations  represent  a  matrix  whose 
coefficients  are  used  to  find  the  20  eigenvalues  and 
eigenvectors  for  use  in  the  EM. 

This  brings  up  a  major  criticism  of  this  method.  The 
matrix  formed  by  the  above  is  not  symmetric r  banded r  and  is 
not  positive  definite.  These  are  not  desireable 
characteristics  of  a  general  matrix.  This  matrix  required 
the  use  of  the  IMSL  routine  EIGRF  (10 :EIGRF-l-EIGRF-2)  to 
find  the  eigenvalues  and  eigenvectors.  This  routine  is 
particularily  difficult  to  use  and  did  not  interface  well 
with  the  VAX  11/780  computer. 


STUDY 


ARTICLE 


RATIO 


EIGENVALUE 

-0.1399 

-0.1400 

0.9999 

COEFFICIENT 

-821.28 

-832.04 

0.8812 

V(l,l) 

0.0269 

0.0267 

1.0075 

V(l,2) 

0.0520 

0.1516 

1.0078 

V{1,3) 

0.0734 

0.0729 

1.0069 

v(l,4) 

0.0898 

0.0891 

1.0090 

v(l,5) 

0.1000 

0.0993 

1.0070 

v(l,6) 

0.1034 

0.1027 

1.0068 

v(l,7) 

0.1003 

0.0996 

1.0070 

v(l,8) 

0.1416 

0.1406 

1.0071 

v(l,9) 

0.1731 

0.1719 

1.0070 

v(l,10) 

0.1927 

0.1913 

1.0073 

v(l,ll) 

0.1994 

0.1979 

1.0076 

v(l,12) 

0.1998 

0.1984 

1.0075 

v(l,13) 

0.2440 

0.2423 

1.0070 

v(l,14) 

0.2714 

0.2695 

1.0071 

v(l,15) 

0.2806 

0.2786 

1.0072 

v(l,16) 

0.2976 

0.2955 

1.0071 

v(l,17) 

0.3304 

0.3281 

1.0070 

v(l,18) 

0.3410 

0.3386 

1.0071 

v(l,19) 

0.3654 

0.3628 

1.0073 

v(l,20) 

0.3747 

0.3720 

1.0073 

Table  I  lists  a  partial  conparison  of  the  output  of  the 
IMSL  routine  EIGRF.  This  is  only  a  sample  of  the  complete 
output  and  clearly  shows  the  results  (19:414)  Shih  and 
Skladany  achieved  are  correct.  The  eigenvectors  have  an 
average  error  -  0.728%.  The  eigenvalues  were  the  same  as 
the  article  (within  4  decimal  places) .  The  coefficient  is  a 
little  different  from  the  Shih  and  S)cladany  coefficient 
because  of  the  slight  differences  in  the  system  of 
eigenvectors.  Despite  the  difference,  the  errors  are 
compensating  and  account  for  a  negligible  difference  in  the 
temperature  histories.  The  coefficients  were  calculated 
using  the  initial  condition  and  solving  the  resulting  matrix 
by  the  IMSL  routine  LEQT2P  (10:LEQT2F-l-LEQT2P-3;  24).  All 
IMSL  routines  used  in  this  study  can  be  found  in  Appendix  C. 

RECURSION  FORMULAE 

With  the  eigenvalues,  eigenvectors,  and  multiplying 
coefficients  known,  the  EM  can  be  used  to  solve  the  problem 
(19:413).  The  nodal  point  temperature  histories  can  be 
found  by 

n 

=  X  +  r  (36) 

Table  II  presents  a  partial  comparison  of  the  results 
of  Eg.  (36) .  Temperatures  at  point  Tg  are  used  as  a  sample 
of  the  entire  grid. 


TBAR 

STUDY 

EXACT 

RATIO 

1.0 

35.211 

32.248 

1.092 

2.0 

70.714 

69.683 

1.015 

3.0 

95.422 

94.984 

1.005 

4.0 

113.134 

112.830 

1.003 

5.0 

126.542 

126.216 

1.003 

6.0 

137.181 

136.786 

1.003 

7.0 

145.913 

145.447 

1.003 

00 

. 

o 

153.244 

152.720 

1.003 

9.0 

159.485 

158.919 

1.003 

10.0 

164.845 

164.254 

1.004 

The  above  table  shows  the  comparison  at  point  with 

r*'  -  1  as  compared  to  the  exact  solution  at  the  same  node 
point.  When  plotted  as  a  graph,  these  values  appear  in 
Figure  3  on  the  following  page.  From  Table  II  and  Fig.  3, 
it  is  seen  that  the  eigenvalue  method  compares  very  well 
with  the  exact  solution  after  t  of  2  or  more.  After  this 
time  the  eigenvalue  method  was  within  ±0,3%  of  the  exact 
answer.  Before  t  »  2  the  method  differed  considerably. 
Shih  and  Skladany  offer  an  explanation  for  this  (19:417). 


RATURE 


1 


Figure  3. 
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"Since  initially  the  heating  rate  OT/K.  and  the 
change  of  the  heating  rate  are  large ,  and 

since  ^At  is  related  to  the  space  derivatives  by 

drrAi  =  oi  v*T 

it  follows  that  initially  the  value  of  must  be 
large  because 

Now  we  already  know  that  the  fourth  derivative  is 
embedded  in  the  leading  truncated  term  in  the 
second  order  accurate  finite  difference 

approximation . " 

This  explanation  also  implies  the  truncation  error  will 
be  large  in  the  initial  heating  in  all  three  schemes. 

Equation  (36)  was  used  as  the  recursion  formula  for  the 
EM.  Using  the  notation  of  Figxire  1,  the  following  recursion 
formulae  were  used  in  the  study  as  the  comparison  routines. 
The  explicit  method  used 


Ip-  ((xAi/f?)LTs+Tw4T^T,?-<)Tp*3+  (37) 


The  implicit  method  used 

-^4  2\  I — -I 

Tp  “  -(o(Ai/h)Cls  +^5 


(38) 


and  the  Crank-Nicolson  used 


if'  =  [c  i-2r) /am)lTp  +  fCV’/2n+2cQ 

Qe +T,J“+TN’^+T?4TEViv+Ti*+TjJ?  09) 

where 

h  =r  cs^'A't  /yf' 


ERROR  EVALUATIONS 

An  article  written  by  Towler  and  Yang  was  the  basis  for 
selection  of  an  error  routine  (23:1021-1024).  The  desire 
for  some  overall  way  of  determining  a  method's  accuracy  was 
used  as  a  criteria  in  the  selection  of  the  Root  Mean  Squared 
(RMS)  method.  Stated  simply ,  this  finds  the  difference 
between  the  temperature  at  each  node  and  the  exact 
temperature,  squares  each  point's  differences,  adds  up  all 
of  these,  divides  by  the  total  number  of  inner  node  points, 
and,  finally,  talces  the  square  root  of  the  answer.  Written 
mathematically  as  (23:1023) 


This  provides  a  benchmark  for  each  routine  to  base 
accuracy  evaluations  upon. 

The  exact  solution  found  In  Appendix  A  Is  In  the  form 
of  a  summation  of  terms  Involving  exponentials ,  sines,  and 
cosines.  When  using  normalized  time  (t  =  1,  2,  3,  ...)  the 
exponential  term  will  drive  the  magnitude  of  the  Individual 
term's  contribution.  Thus,  when  n,  m,  and  o  become  large, 
the  term  will  be  small  and  add  little  to  the  summation. 

When  using  normalized  time  it  is  found  that  only  15  terms  In 
the  sum  are  required  for  convergence  to  10"'®  accuracy. 

This  degree  of  accuracy  Is  seldom  necessary  In  light  of  the 
accuracy  of  such  measured  things  as  thermal  conductivity, 
specific  heat,  etc.  To  retain  a  reasonable  amount  of 
accuracy  of  3  or  4  significant  digits,  only  11  terms  are 
required  In  the  summation.  This  translates  Into  significant 
time  savings  In  the  computation  of  the  exact  solution  for 
comparison  purposes. 

RESULTS  OF  2^^  PROBLEM 

The  sample  problem  in  Shih  and  Skladany’s  article  was 
symmetric  in  the  x-y  plane  (19:413).  For  this  reason,  they 
elected  to  take  full  advantage  of  symmetry  and,  thus,  make 
the  computational  domain  fairly  small  (a  20X20  matrix) . 

Full  symmetry  In  everyday  engineering  problems  Is  Idealistic 
(22:127).  Without  symmetry,  the  domain  will  be  analyzed  by 
a  full  grid  with  pitch  h.  This  study  ran  the  test  problem 
both  ways.  Symmetry  was  used  for  the  Initial  verification. 


The  same  problem  was  then  solved  using  the  full  grid  (same 
h)  without  taking  advantage  of  the  symmmetry.  This  approach 
brought  out  some  interesting  observations. 


Taking  advantage  of  synunetry  produces  a  matrix  that  is 
non-symmetric /  non-banded,  and  is  non-positive  definite. 

This  type  matrix  required  the  use  of  a  more  time  consuming 
routine  (EIGRF) .  It  also  required  more  core  memory  as  the 
matrix  must  be  stored  in  its  full  form.  The  full  grid 
analysis  produces  a  banded,  symmetric,  non-positive  definite 
matrix.  This  saved  central  processor  time  (CPU  time)  and 
core  memory  as  only  the  lower  co-diagonals  needed  to  be 
stored  to  compute  the  eigenvalues/vectors. 

Shih  and  Skladany  did  not  do  a  time  comparison  to 
determine  if  the  EM  approach  was  faster  (and  therefore,  less 
expensive)  on  computers  than  the  other  types  of  solving 
schemes.  This  was  accomplished  in  this  study  using  the  VAX 
11/780  and  the  UNIX  command  TIME  (16:473-474).  This  command 
returned  several  pieces  of  information  including  wall  clock 
time  and  the  CPU  time. 


Taking  advantage  of  TIME  and  the  RMS  error  routine,  the 
four  methods  were  run  in  symmetry  and,  then,  in  full  grid 
forms.  The  results  were  then  compared  to  the  exact 
solution. 
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TIME  COMPARISONS 


The  following  tables  represent  the  time  comparisons  of 
the  specific  methods.  The  CPU  times  are  accurate  to  within 
1/10  second  and  the  wall  clock  times  are  accurate  to  within 
1/60  second  (16:473-374).  The  wall  clock  times  may  not  be 
consistent  because  runs  made  during  different  times  of  the 
day  took  more  or  less  wall  clock  time  according  to  the 
computer's  load.  The  CPU  times  will  be  consistent  as  this 
Is  the  time  the  central  pocessor  takes  to  execute  the 
program . 


TABLE  III 

Time  Comparison  of  Methods  In  Symmetry. 


METHOD 

CPU (SEC) 

WALL (HR: MIN: SEC) 

EIGENVALUE 

0.8 

0:00:42 

EXPLICIT 

0.7 

0:02:13 

IMPLICIT 

4.2 

0:16:34 

CN 

6.4 

0:23:56 

Here  It  Is  easy  to  see  the  EM  Is  very  comparable  to  the 
explicit  method  and  Is  clearly  faster  than  the  last  two 
methods.  All  programs  were  run  to  achieve  similar 
accuracies  (three  significant  digits  the  same) . 


METHOD 

CPU (SEC) 

WALL(HR:MIN:SEC) 

EIGENVALUE 

40.1 

3:04:19 

EXPLICIT 

15.1 

1:44:47 

IMPLICIT 

6.1 

0:41:41 

CN 

26.4 

1:29:00 

The  value  of  the  use  of  the  full  grid  approach  becomes 
apparent.  When  comparing  Tables  III  and  IV  the  CPU  times  do 
not  seem  consistent.  The  explicit  time  in  Table  IV  can  be 
explained  by  realizing  the  time  restriction  in  Eg.  (7)  must 
be  satisfied.  In  this  case  At  >=  1/1000  to  achieve  the 
necessary  accuracy.  This  requires  1000  Incremental 
calculations  to  get  to  t  =  1  I  The  implicit  time  is  less 
than  explicit  because  it  does  not  require  this  sort  of  time 
step.  To  achieve  the  same  accuracy  the  At  1/50  .  The  CN 
scheme  takes  more  time  in  both  tables  and  is  therefore 
consistent. 

The  problem  of  having  to  find  the  eigenvalues  and 
eigenvectors  starts  to  become  clear  here.  Finding 
eigenvalues /vectors  in  full  grid  results  in  nxn  matricles 
being  formed.  The  full  grid  sample  used  121  Internal  mesh 
points  (11X11  grid,  121X121  matrlxl).  In  full  grid  mode  the 
method  now  takes  4  times  as  long  to  run.  In  this  case  all 


of  the  other  methods  are  preferable  to  the  EM;  just 
considering  the  time  required  to  execute. 


-  EIGEN 

0 — 0  EXPLICIT 
B — B  IMPLICIT 
^ - H  CN 


EXTENSION  OP  SAMPLE  PROBLEM  INTO  3-D 


A  literature  search  of  the  numerical  techniques  in 
solving  the  heat  equation  revealed  very  few  sources  that  go 
beyond  a  two  dimensional  case  (17:45-48;  4:142-144; 
3:162-171).  It  iS/  therefore,  instructive  to  examine  a 
problem  in  three  dimensions. 

Several  questions  arrise  immediately  when  considering 
the  three  dimensional  case. 

1.  Does  an  analytical  solution  exist? 

2.  What  boundary  conditions  are  appropriate? 

3.  Does  the  method  accuracy  change  in  3-D? 

4.  How  much  more  time  is  required  to  run? 

Answers  to  these  questions  are  supplied  in  the 
following  sections. 

THE  THREE  DIMENSION  PROBLEM  STATEMENT 

The  simplest  case  to  study  is  a  cube  of  dimension  L  on 
a  side.  A  more  general  problem  also  worth  considering  is  a 
parallelpiped  (a  box)  with  dimensions  as  show  in  Figure  4. 

The  general  case  can  be  easily  modified  by  ma)cing  the 


lengths  a,  b,  and  c  all  equal  to  L  and  the  cube  case  is 
back.  The  program  written  took  advantage  of  this,  thus 
saving  work  and  time. 


FIGURE  4.  Three  dimension  test  problem 


All  the  surfaces  of  the  box  will  be  at  the  same  ■ 
200  All  the  interior  points  are  at  0  initially. 

Therefore,  the  governing  equation,  boundary  conditions,  and 
initial  condition  are 

=  k  («) 

Where 

T(0,y,2,t)  «  To 
T(a,y,z,t)  »  To 
T(x,0,z,t)  «  To 
T(x,b,z,t)  «  To 
T(x,y,0,t)  -  To 
T(x,y,c,t)  ■  To 


T(x,y,z,0)  -  T;^ 

Using  the  same  normalizations  as  the  2-D  case  and 
refer ing  to  Appendix  A,  the  solution  to  Eg.  (41)  is  given  by 
Eg.  (A-41) .  Indeed,  the  boundary  conditions  are  appropriate 
and  there  is  an  analytical  solution  to  compare  the  methods 
against. 

Time  restrictions  during  the  study  only  allowed  for  the 
implicit  method  to  be  developed  in  three  dimensions  as  a 
comparison  routine  to  the  exact  and  eigenvalue  methods. 
Little  is  lost  here  because  the  2-0  results  extend  into  the 
3-D  case. 

Recalling  that  the  number  of  internal  mesh  points 
determines  the  n  dimension  of  the  eigenvector/value  solving 
routine,  a  balance  must  be  met.  Enough  mesh  points  must  be 
supplied  to  Insure  accuracy  (not  necessarilly  stability) 

(15)  while  at  the  same  time  too  many  mesh  points  cause  the 
problem  to  run  extremely  long.  In  the  eigenvalue  case,  a 
cube  with  8  divisions  (7  internal  points)  on  a  side  was  the 
maximum  number  the  VAX  computer  could  handle  before  running 
out  of  core  memory.  The  core  memory  reguires  bytes  (10) 
available  for  the  eigenvalue  method  to  use  even  in  full  grid 
form  (banded  storage  mode) . 

The  EM  is  not  an  ideal  method  if  computer  memory  space 
is  at  a  premium.  This  is  another  criticism  of  this  method. 
Shih  and  Skladany  did  not  explore  past  the  20X20  matrix  they 
used  and  this  aspect  of  the  eigenvalue  method  was  not  made 
apparent.  This  could  present  a  definite  problem  to  the  user 


if  he  is  trying  to  solve  a  heat  problem  for  something  large 
such  as  a  nuclear  reactor  core. 


Figures  5  and  6  are  the  plots  of  the  EM  and  IM  CPU  time 
requirements f  respectively.  It  appears  the  EM  is 
exponential  in  nature  with  respect  to  the  number  of  mesh 
points  while  the  IM  is  simply  linear.  The  curves  drawn  in 
the  figure  represent  the  mean  value  of  the  data  points. 

The  first  curve  (Fig.  5)  was  dravm  by  developing  the 
following  data: 

TABLE  VII 

Data  Used  for  the  Determination  of  Kl. 


GRID 

1 

POINTS 

TIME (SEC) 

3X3X3 

27 

5X4X3 

60 

5.0 

4X4X4 

64 

6.2 

11X11 

121 

40.1 

5X5X5 

125 

55.4 

The  curve  in  Fig.  5  is  a  plot  of 


TIME  s  exp  (KlxN) 


where 


Kl  «  an  averaged  constant 
n  «  number  of  mesh  points 


The  value  of  K1  was  arrived  at  by  a  simple /  although 
somewhat  crude r  technique  of  averages.  Five  runs  were  made 
and  timed,  each  with  different  numbers  of  mesh  points.  In 
each  case  the  log  of  the  time  was  divided  by  the  number  of 
mesh  points.  These  five  numbers  are  then  averaged  and  yield 

K1  -  0.0317 

The  behavior  of  the  implicit  with  respect  to  time 
required  appears  linear.  Assuming  this  to  be  true,  the 
governing  equation  is  simply 

TIME  «  K2xN  +  I  (43) 

where 

I  *  y  intercept 

A  linear  regression  produced  a  slope  of  K2  *  0.0472 
with  an  Intercept  of  I  =  0.4555  (18) .  The  data  used  in  the 
plot  are 


GRID 

POINTS 

3X3X3 

27 

5X4X3 

60 

4X4X4 

64 

11X11 

121 

5X5X5 

125 

TIME (SEC) 


The  exponential  time  requirement  is  a  definite  handicap 
for  the  EM  unless  a  coarse  mesh  is  used  in  the  original 
problem. 


ACCURACY  OF  3^0  RESULTS 

Little  was  learned  concerning  the  accuracy  of  the 
methods  except  that  the  overall  accuracy  of  both  methods 
improved  with  the  addition  of  more  Internal  mesh  points. 

This  is  not  an  unexpected  result,  but  gratifying  to  know  the 
theory  correctly  predicts  this.  Table  IX  lists  the  results 
of  the  accuracy  run  for  the  7X7X7  grid  cube. 


METHOD 

TBAR 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

EM 

2.32 

0.15 

0.11 

0.03 

0.01 

0.00 

0.00 

0.00 

0.00 

0.00 

IM 

3.63 

1.79 

0.73 

0.21 

0.03 

0.10 

0.12 

0.11 

0.09 

0.07 

Again,  as  in  the  2-D  case,  the  EM  is  more  accurate  that 
the  IM.  The  accuracy  seems  to  get  better  after  TEAR  greater 
than  5.0  .  The  problem  now  has  temperatures  on  all  surfaces 
which,  in  turn,  causes  a  rapid  heating  of  the  interior.  The 
material  will  therefore  reach  a  steady  state  condition 
faster  than  in  the  two  dimensional  case.  This  would  account 
for  the  apparent  exactness  of  the  answers  at  late  times. 

NON-SYMMETRIC  3-D  PROBLEM 

An  additional  3-D  problem  was  run  to  deliberately 
remove  any  gain  from  symmetry.  Very  little  was  gained  from 
the  exercise,  except  the  knowledge  the  method  still  works. 

The  problem  was  to  make  a  parallelpiped  having 
dimensions  such  that 

a  =  5 
b  =  4 
c  =  3 


This  box  shape  would  make  it  very  hard  to  take  advantage  of 


symmetry.  Since  the  other  problems  were  run  with  full  grid 
mode,  this  was  a  trivial  problem  to  run. 

The  accuracy  results  were  as  follows 


TABLE  X 

Comparison  of  Accuracy  for  5X4X3  Parallelpiped. 


METHOD 

1 

2 

3 

4 

5 

TEAR 

6 

7 

8 

9 

10 

EM 

1.43 

0.07 

0.07 

0.04 

0.02 

0.01 

0.00 

0.00 

0.00 

0.00 

IM 

1.15 

0.14 

0.16 

0.08 

0.03 

0.01 

0.00 

0.00 

0.00 

0.00 

The  same  behavior  is  exhibited  here  as  in  the  cube  in 
Table  IX.  The  answers  do  not  appear  symmetric  because  of 
the  lack  of  symmetry  in  the  region,  as  expected.  The 
answers  also  appear  to  be  more  accurate  in  Table  X  than 
Table  IX.  This  is  due  in  part  to  less  numbers  of  node 
points.  The  3-D  geometry  caused  a  good  increase  in  accuracy 
over  the  2-D  case,  but  the  accuracy  should  improve  a  little 
if  not  as  many  calculations  (roundoff  error)  are  made. 


V.  EVALUATION  OF  THE  EIGENVALUE  METHOD 


Based  upon  results  of  the  previous  chapter,  an  overall 
rating  is  made  of  the  eigenvalue  method. 


THE  STRONG  POINTS 

The  accuracy  of  the  EM  is  better  than  the  other  methods 
(EX,IM,CN) .  Even  at  earlier  times  the  EM  is  still  strong. 
This  better  accuracy  remains  throughout  a  transient, 
regardless  of  duration. 

The  reason  given  for  the  improved  accuracy  is  the 
exactness  of  the  time  derivative  in  the  heat  equation.  The 
other  methods  are  similarily  derived  with  some  form  of 
discretized  approximation  via  Taylor  series  (2:347-348). 

The  EM  is  even  more  accurate  then  the  best  of  the  comparison 
methods,  CM,  in  the  time  derivative. 

The  appearance  that  the  other  methods  approach  the  same 
accuracy  of  the  EM  at  later  times  is  misleading.  Figure  4 
shows  the  error  of  the  comparison  methods  to  converge  to 
that  of  the  EM.  Mathematically,  we  know  this  to  be  untrue. 
The  difference  between  the  EM  and  the  other  methods  accuracy 
should  remain  constant.  This  convergence  is  brought  about 
by  the  region  under  consideration  approaching  the  steady 
state  condition. 

The  efficiency  of  the  EM  must  lie  on  the  side  of  being 


an  asset.  Consider  a  hypothetical  problem.  A  small  nuclear 
reactor  (as  on  a  submarine)  has  a  problem  requiring  a 


shutdown  and  subsequent  cooling.  Technical  problens  dictate 
a  slow  cool  down  taking  1  week.  What  temperature  behavior 
happens  1  meter  from  center? 

Using  the  EM  and  assuming  a  1-0  treatment  (the  reactor 
is  a  symmetric  cylinder  with  variation  only  in  the  radial 
direction) ,  how  long  would  it  take  the  computer  to  determine 
reactor  temperatures  at  the  point  at  t  «  1  week? 

Assume  each  method  divides  the  1  meter  length  into  100 
intervals  for  accuracy  and  allows  the  reference  time,  t^,  to 
be  1  second.  Using  Fig.  (5) «  we  find  the  EM  requires 
approximately  20-25  seconds  of  CPU  time  to  solve  the 
problem.  This  time  includes  determination  of  the 
eigenvalues/  eigenvectors/  and  coefficients.  Figure  6  shows 
that  it  takes  the  IM  approximately  5  seconds  of  CPU  time  (5 
sec  to  reach  T  •  10) .  This  means  it  requires  1/2  second  to 
reach  t  «  1.0.  Since  the  reference  time  is  unity/  T  «  t 
and  the  IM  must  Iterate  till  it  reaches  1  week/  or  604/800 
seconds.  The  IM  would  then  require  302/400  CPU  seconds  to 
reach  the  same  solution  found  by  the  EM  in  20-25  seconds  CPU 
time.  There  is  little  doubt  which  scheme  the  engineer  on 
the  submarine  will  choose. 

Two  more  assets  are  possessed  by  the  eigenvalue  method. 
The  first  is  the  unconditional  stability  of  the  method 
brought  about  because  of  the  analytical  nature  in  the 
temporal  domain.  The  second  point  is  the  omission  of  most 
of  the  terms  in  Eq.  (36).  This  is  allowed  because  the 
coefficients  become  so  small  in  some  cases  that  most  terms 


contribute  very  little  to  the  overall  answer  and  can  be 
eliminated  without  affecting  the  answer.  Undoubtedly/  a 
solution  requiring  5  terms  runs  faster  than  one  requiring 
100.  This  point  is  moot.  The  time  saved  in  eliminating  the 
extra  terms  is  often,  as  in  the  case  of  large  numbers  of 
mesh  points,  negligible  when  compared  to  the  time  required 
to  find  all  the  eigenvalues,  eigenvectors,  and  coefficients. 

THE  WEAK  POINTS 

ParadoxicaMy,  one  of  the  EM's  strongest  points  can 
also  be  its  weakest  point.  It  has  been  demonstrated  the  EM 
performs  poorly  with  large  numbers  of  mesh  points  (requires 
large  CPU  times) .  This  may  not  be  an  Insurmountable  problem 
if  one  is  using  a  very  capable  computer,  but  when  restricted 
to  a  machine  such  as  the  VAX  11/780,  it  is  a  definite 
problem. 

The  heart  of  the  problem  is  the  requirement  to  find 
complete  sets  of  eigenvalues,  eigenvectors,  and 
coefficients.  This  study  used  standard  IMSL  routines. 

These  routines  are  widely  available  and  are  commonly  used  in 
such  studies  (20;  23) .  The  EM  has  potential  of  being  very 
powerful  if  better,  faster  routines  were  made  available  to 
find  these  eigenvalues /vectors.  An  extensive  literature 
search  (1;  3;  5;  6;  7;  11;  17) ,  indicates  this  problem  is 
being  given  top  priority  in  the  numerical  analysis  community 
(12). 

Another  weak  point  of  the  EM  is  the  inability  to  hand 


check  any  of  the  calculations.  The  nature  of  finding 
eigenvalues  and  eigenvectors  makes  for  extreme  difficulty  in 
a  hand  check  of  the  calculations,  especially  for  a  large 
matrix.  Imagine  trying  to  hand  calculate  a  100X100 
eigenvector /value  problem.  This  discourages  the  novice  user 
from  taking  advantage  of  the  EM.  If  the  problem  has  no 
analytical  solution  or  is  not  easily  approximated,  the  EM 
user  will  have  little  confidence  in  his  answers. 


VI.  CONCLDSIONS  AND  RECOMMENDATIONS 


CONCLUSIONS 

The  eigenvalue  method  is  a  good  numerical  technique 
worthy  of  inclusion  in  any  problem  solver's  collection  of 
numerical  techniques.  It  is  not  a  panacea  capable  of 
solving  problems  other  methods  such  as  explicit,  implicit, 
or  Crank-Nlcolson  cannot. 

Like  all  things,  the  EM  has  its  good  points  and  its 
bad  points.  It  should  be  considered  better  in  accuracy  than 
most  discretization  methods.  It  is  expecially  well  suited 
in  large  parametric  studies.  It  does  not  lend  itself  well 
to  large  spatial  studies  involving  great  numbers  of  mesh 
points.  The  improved  accuracy  and  the  slower  execution 
times  can  be  paired  against  one  another  to  achieve  a 
reasonably  fast  and  accurate  method. 

When  used  as  a  general  problem  solver  the  EM  is  about 
the  same  as  any  other  technique.  It  is  when  used  in  large 
transient  time  applications  that  the  EM  is  vastly  superior. 

The  EM  is  unconditionally  stable  because  of  the 
exactness  of  the  time  derivative  of  the  temperature  in  the 
heat  equation. 

The  EM  translates  easily  into  three  dimensions.  It 
follows  the  method  should  also  work  well  in  other  coordinate 
systems,  such  as  cylindrical  or  spherical. 

There  is  a  potential  of  high  CPU  useage  time  in  the  EM 


if  large  numbers  of  internal  mesh  points  are  encountered. 


The  EM  will  always  be  slightly  better  in  accuracy  than 
the  standard  discretization  techniques  used  in  this  study. 
This  accuracy  improvement  is  made  available  by  the 
analytical  nature  of  the  time  derivative. 

RECOMMENDATIONS 

It  is  recommended  that  a  more  detailed  investigation  be 
made  into  the  use  of  the  eigenvalue  method.  Extensions  of 
the  method  in  spherical  and  cylindrical  coordinates  should 
be  made.  Also  of  interest  would  be  how  the  EM  performs  on 
regions  with  irregular  boundaries. 

It  is  further  recommended  that  a  better  description  of 
the  amount  of  CPU  time  required  to  execute  the  programs  be 
given.  The  one  expressed  in  this  study  used  only  five  data 
points.  More  sophisticated  curve  fits  should  be  employed  to 
verify  or  refute  the  exponential  nature  of  the  time 
requirement  of  the  EM. 

Another  recommendation  is  the  use  of  the  EM  as 
presented  in  this  study  with  faster r  more  capable  machines. 
It  would  be  interesting  to  reaccomplish  this  study  on  the 
HARRIS  800 r  CYBERr  or  CRAY  computers. 

It  is  recommended  that  a  better  method  be  determined  to 
find  eigenvalues  and  eigenvectors.  If  such  methods  do  not 
yet  exist,  then  it  is  suggested  that  carefull  review  of  the 
existing  methods  (IMSL)  be  made  with  the  goal  of  improving 
their  speeds  and/or  accuracies.  Also  in  this  area,  it  is 
recommended  that  a  carefull  investigation  be  made  to 


determine  if  all  the  eigenvalues  and  eigenvectors  need  to  be 
calculated  In  order  to  find  the  correct  multiplying 
coefficients  for  the  solution  of  the  problem. 

It  Is  recommended  that  a  study  be  made  to  determine  a 
way  to  predict  the  balance  between  the  EM  accuracy  and  mesh 
point  spacing.  This  would  save  much  time  In  trial  and 
error.  Once  this  balance  has  been  determined  it  would  be 
instructive  to  make  a  computer  study  to  further  explore  long 
transient  problems. 

It  Is  recommended  that  an  investigation  of  the  method 
of  lines  (Runge-Kutta)  be  made  and  results  compared  to  the 
eigenvalue  method  to  identify  similarities  and  differences. 

A  time  comparison  between  the  two  methods  would  prove 
extremely  useful. 

Finally,  it  is  recommended  that  problems  involving 
radiative  surfaces,  heat  sources,  and  differing  material 
properties  (such  as  specific  heat,  thermal  transmissivity, 
etc.)  be  studied  using  the  EM.  Furthermore,  differing 
boundary  conditions  such  as  Neuman,  Dirlchlet,  and  Robins 
should  be  examined  and  evaluated. 


APPENDIX  A 


DERIVATION  OF  THE  ANALYTICAL  SOLUTION 

In  this  appendix,  verification  of  the  sample  problem 
presented  by  Shlh  and  Skladany  will  be  accomplished 
(19:411).  A  three  dimensional  problem  is  also  solved  in 
this  report.  Therefore,  it  is  convenient  to  extend  the 
derivation  to  three  dimensions.  This  was  not  done  in  the 
previously  mentioned  article. 

THE  SAMPLE  PROBLEM 

Given  the  heat  equation 

with  boundary  conditions  and  initial  conditions  as 

T(0,y,z,t)  «  To 
T(a,y,z,t)  »  •So 
T(x,0,z,t)  =  To 
T(x,b,z,t)  =  To 
T(x,y,0,t)  »  To 
T(x,y,c,t)  «  To 
T(x,y,z,0)  »  Ti 

Now,  let  the  normalized  temperature  be  defined  as 


Then  the  normalized  boundary  conditions  are  now 

B  (0,y,z,t)  -  0 
Q  (a,y,z,t)  *  0 
@(x,0,z,t)  =  0 
© (x,b,z,t)  *  0 
&  (x,y,0,t)  =  0 
© (x,y,c,t)  a  0 
0 (x,y,z,0)  -  1 

In  terms  of  the  normalized  temperature,  Eq.  (A-1)  now 
appears  as 

<»-3 

where 

c/-  k/fCv 


SEPARATION  OF  VARIABLES 

Assuming  the  method  of  separation  of  variables 
technique  is  valid  (because  of  the  same  temperature  on  all 
surfaces,  then  there  exists  some  solution  (4:34) 

Y(ij^*ZCz)  (A-4 

and  differentiating  with  respect  to  the  appropriate 
independent  variable 


=  XYZT' 

,afe/9x^  =  x"rzT 


(A-5) 


(A-6) 


(A-7) 


cte/^^=XY2"T 


(A-8) 


Using  these,  Eq.  (A-3)  becomes 


XYZT'  =  o<CXl'  YZT  +  XY*ZT  +  XYZ-’T'D  ,»-9, 


dividing  both  sides  by  XYZT 


Y/t  =  ocC  xVx  +  Y"A+  Z'Vz'l  (»-j 


rearranging  gives  the  form 


r'/kT  -  Y'W-  2"  /2  =  XVX 


(A-11 


In  typical  separation  of  variables  poblems,  both  sides 
must  also  equal  some  constant,  say, 


xvx  =  -A 


(A-12 


considering  the  x  dependence 


This  equation  typically  has  solutions  of  the  form 


X(iO  =  Cl  5in  vTX  +  Ce  eos  \/J  X 

applying  the  normalized  boundary  conditions 

Yet)  -  Ci  s\r\  M  +  Ca  Cos^oi 

for  this  to  be  true  C2.  =  0.  AlsOf 

Xca1  =  CiSinN/7a.  =  c> 

to  avoid  the  trlval  solution  ^  0,  therefore^ 

\/Ja  =  nfl 
h  -  (  nn  /a) 


with 


l!  Cft  Sin  (rnrx/a) 


where 


similar  treatments  of  the  y  and  z  variables  will  yield 


^in  CmTTy/b^ 

Zj  Q)  srn  (.ott'2 

0=1 

where  n,  m,  and  o  are  odd  as  from  before.  The  time 
dependence  Is  found  from 

k(  TT”  —  o 

X  separation  constant 
y  separation  constant 
z  separation  constant 

k  (  y6rf\4  Yft  ^  fr\f*to 

With  these  sxibstitutlons ,  Eq.  (A-22)  becomes 

-p  J^Artte  1  ^ 

whose  solution  has  the  form 

Tl-b)  =  3  ®(pC-  frvo.-t) 


T'4 

where 

An  = 
/Son  * 

and  let 


applying  the  Initial  condition 


Tco^  =  <i  exp  co^  -  1 


(A-25 


which  implies  that  3=  1.  Therefore, 


"110  =  exp 


(A-26 


Combining  all  the  constants  thus  far  as 


and  using  Eqs.  (A-19) ,  (A-20) ,  (A-21) ,  and  (A>24)  and 
recalling  Eq.  (A-*3)  one  gets  the  normalized  temperature 
solution 


00  0000 

Xnmo  expG* 

'  rvim»io»i  ~ 


4-  /b^  +  oVc^)  kll  *  -il  nc  hTIV/a) 

Sih(mTrvi/b^  iiniotrz/c'i 


FOURIER  SERIES 

Let  f(Xry,z,0)  >  1  .  From  Fourier  series,  any 


function  may  be  expanded  as  a  sine  series  (4:110) 


(A-29) 


QO 


2  An  sin  (nnx/cii 


and 


a 

An  =  (■’Zt'a)  S,-fcx,  Li,z>  -?in  ( nirx/o^  dv 
An  =  (z/q)£  1  5in  CnTTvA^clx 


(A-30) 


(A-31) 


but,  A^  can  also  be  expanded  as 


QO 


At,-  Z  Bnn,  S(n(mTry/b^ 


where 


I 


(A-32) 


Bnm  =  J  An^in(  m¥3/b)oly  '*-33' 


and,  like  above,  can  also  be  expanded  as 

oo 

■'nm  “ 


Snm  ^  Knmo  'Sin  (oTrz./c'\ 


{A-34) 


0=1 


where 


Xnmo  =  (a/c)  I  B  ^lacoTTz/c^  6z 


(A-35) 


plugging  Eqs.  (A-31)  and  (A-33)  into  Eq.  (A-35)  yields 


57 


A  mo 


abc 

=  (S/abe)  nf  sinCnnx/a)  sin^m/ryA^ 

Sir\(OTrz/c)clzcl}jdi(  ,*-36) 

The  integral 

L 

J  'Sin  (  (A-37) 

has  value  only  when  n  is  odd  and  is  equal  to 

2L  /  mr  (A-38) 

With  this  integral  appearing  three  times  in  Eq.  (A-36)  the 
constant  is  clearly 

KAmo  “  (  (a-39) 

Homo  ~  /  fnmoTT^'^ 


THE  NORMALIZED  SOLDTION 

Once  the  value  for  the  constant  is  kno%nir  it  Is  used  in 
Eq.  (A-28)  to  give  the  result 

-  (rf/Q^+ 

yinfirMiy /Jilsinfu'ii y/b^SmtoHi/c) 

This  is  an  extended  version  of  the  solution  found  in  the 
Shih  and  Skladany  article  and  verifies  their  exact  solution 
(19:410).  This  solution  is  valid  only  for  terms  where  n,  m, 
and  o  are  odd  and  not  equal  to  zero. 


APPENDIX  B 


APPLICATION  OF  THE  EIGENVALUE  METHOD 
TO  A  SMALL  SYSTEM  OP  EQUATIONS 


The  core  of  the  eigenvalue  method  is  the  vector 
representation  of  the  heat  equation  as 


where 


’  du 

diz 

dsi 

dit 

(B-2) 


Shlh  and  Skladany  applied  the  above  to  a  2x2  matrix 
system  (19:412).  This  appendix  will  do  the  same/ 
paralleling  their  development  with  a  few  extra  steps  allowed 
for  clarity. 

Temporarily  dropping  the  vector  and  matrix  notation  and 
looking  only  at  the  equations  of  the  sample  matrix,  Eq. 

(B-1)  can  be  written 


0(t)  -  D  =  O  <’ 

this  type  differential  equation  is  recognized  as  having 
solutions  of  the  form 


=  XieicpOit)  +  yiOcpCAtO 


(B-4) 


=  )4ej(p(Ait)+  (B-5) 

differentiation  with  respect  to  time 

=  >(X.expO,-tU  AsL4i®9<'Aat')  ‘®"'’ 

0£rt)  =  ^iXtexprMH  A^Jj^eKpfA^i) 

plugging  Eqs.  (B-6)  and  (B-7)  into  (B-4)  and  (B-5)  gets 

XiXiexp(^tt)+  X2yie/p(^jt)=  d„[(X,expCXft) 

+ ckK  Xz-tJ]  4-c)ieC(Xj,e<p(  A(tl+yjeK|>Oit^  '®-“> 
AXtecpCAttlf  AzyeescpCXst) »  ctiC(fX,exp(A(i) 

4- u,e)cpOztl^+d2  OXieKpCkU  y£exj)(Xrl)3  <®-®> 

placing  everything  on  the  right  side  and  simplifying 

>- exp(A,t^(diiXi4dft4-AiXiHexp(AztXdii\ji+tl«yt-Aiyi) 

1=  sxpCXii^iVi-kla^- 


Equations  (B-10)  and  (B-11)  when  written  in  matrix  form 


0=  eKpa't)[<4,  cfe->,JUJ+afl.ait)|d„ 


from  Eq.  (B-2) 


which  implies  li  and  It  are  the  eigenvalues  of  and  Cifi  ktl'*’ 
and  C  iji  are  the  corresponding  eigenvectors .  The  terms 

exp(  X|t)  and  exp(>ait)  are  just  multiplying  scalars.  If  Ai  ,  At 
•  CXi  /  and  'jt3^  be  determined,  Eq.  (B-1)  is 


solved. 


APPENDIX  C 


STANDARD  IMSL  ROUTINES  USED  IN  THE  STUDY 


This  appendix  contains  copies  of  the  standard  IMSL 
routines  used  in  the  various  programs  in  the  study.  This 
appendix  is  provided  to  the  reader  to  aide  in  the  evaluation 
of  the  solving  routines  used.  Any  questions  concerning 
method  of  solution  once  the  programs  called  these  routines 
should  be  refered  to  the  subroutine  author. 


IMSL  BOUTINZ  MAME 


EIGRF 


-  EIGENVALUES  AND  (OPTIONALLY)  EIGENVECTORS  OF 
A  REAL  GENERAL  MATRIX  IN  FULL  STORAGE  MODE 

-  CALL  EIGRF  (A,N, lA, IJOB,W, Z , IZ , WK, lER) 

-  THE  INPUT  REAL  GENERAL  MATRIX  OF  ORDER  N 
WHOSE  EIGENVALUES  AND  EIGENVECTORS  ARE 
TO  BE  COMPUTED.  INPUT  A  IS  DESTROYED  IF 
IJOB  IS  EQUAL  TO  0  OR  1. 

-  THE  INPUT  ORDER  OP  TEE  MATRIX  A. 

-  THE  INPUT  ROW  DIMENSION  OF  MATRIX  A  EX.ACTLY 
AS  SPECIFIED  IN  THE  DIMENSION  STATEMENT  IN 
THE  CALLING  PROGRAM. 

-  THE  INPUT  OPTION  PARAMETER.  WHEN 
IJOB  -vO,  COMPUTE  EIGENVALUES  ONLY 
IJOB  -  1,  COMPUTE  EIGENVALUES  AND  EIGEN¬ 
VECTORS  . 

IJOB  »  2,  COMPUTE  EIGENVALUES,  EIGENVECTORS 
AND  PERFORMANCE  INDEX. 

IJOB  ■  3,  COMPUTE  PERFORMANCE  INDEX  ONLY. 

IP  THE  PERFORMANCE  INDEX  IS  COMPUTED,  IT  IS 
RETURNED  IN  WK(1) .  THE  ROUTINES  HAVE 
PERFORMED  (WELL,  SATISFACTORILY,  POORLY)  IF 
WK(I)  IS  (LESS  THAN  1,  BETWEEN  1  AND  100, 
GREATER  THAN  100) . 

W  -  THE  OUTPUT  COMPLEX  VECTOR  OP  LENGTH  N, 

CONTAINING  THE  EIGENVALUES  OF  A. 

NOTE  -  THE  ROUTINE  TREATS  W  AS  A  REAL  VECTOR 
OP  LENGTH  2*N.  AN  APPROPRIATE  EQUIVALENCE 
STATEMENT  MAY  BE  REQUIRED.  SEE  DOCUMENT 
EXAMPLE. 

Z  .  -  THE  OUTPUT  N  BY  N  COMPLEX  MATRIX  CONTAINING 

THE  EIGENVECTORS  OF  A. 

THE  EIGENVECTOR  IN  COLUMN  J  OF  Z  CORRES¬ 
PONDS  TO  THE  EIGENVALUE  W(J) . 

IF  IJOB  «  0,  Z  IS  NOT  USED. 

NOTE  -  THE  ROUTINE  TREATS  2  A  REAL  VECTOR 
OP  LENGTH  2*N*N.  AN  APPROPRIATE  EQUIVALENCE 
STATEMENT  MAY  BE  REQUIRED.  SEE  DOCUMENT 
EXAMPLE. 

12  -  the  input  row  DIMENSION  OP  MATRIX  Z  EX.\CTLY 

AS  SPECIFIED  IN  THE  DIMENSION  STATEMENT  IN 
THE  CALLING  PRCGR.2LM.  IZ  MUST  BE  GREATER 
THAN  OR  EQUAL  TO  N  IF  IJOB  IS  NOT  EQUAL  TO 
ZERO. 

NK  -  WORK  AREA,  THE  LENGTH  OF  WK  DEPENDS 


ON  THE 
IJOB  « 

VALUE  OF  IJOB, 
0,  THE  LENGTH 

WEEN 
OF  WK 

IS 

AT 

LEAST 

N. 

IJOB  « 

1, 

THE 

LENGTH 

OF 

WK 

IS 

AT 

LEAST 

2N. 

IJOB  « 

2, 

THE 

LENGTH 

OF 

WK 

IS 

AT 

LEAST 

(2-l-N)N. 
IJOB  «  3, 

THE 

LENGTH 

OF 

WK 

IS 

AT 

LEAST 

1. 

PURPOSE 

USAGE 

arguments  a  . 

N 

lA 

IJOB 
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SIGRP-l 


lER  -  ERROR  PARAMETER.  (OUTPUT) 

TERMINAL  ERROR 

lER  »  128+J,  INDICATES  THAT  SQRH3F  FAILED 
TO  CONVERGE  ON  EIGENVALUE  J.  EIGENVALUES 
J+l,J+2, . . . ,N  HAVE  BEEN  COMPUTED  CORRECTLY. 
EIGENVALUES  ARE  SET  TO  ZERO. 

IP  IJOB  »  1  OR  2  EIGENVECTORS  ARE  SET  TO 
•  ZERO.  THE  PERFORMANCE  INDEX  IS  SET  TO  1000. 

WARNING  ERROR  (WITH  FIX) 

lER  -  66,  INDICATES  IJOB  IS  LESS  THAN  0  OR 
UOB  IS  GREATER  THAN  3.  IJOB  SET  TO  1. 
lER  ■  67,  INDICATES  IJOB  IS  NOT  EQUAL  TO 
ZERO,  AND  IZ  IS  LESS  THAN  THE  ORDER  OF 

MATRIX  A.  IJOB  IS  SET  TO  ZERO. 

« 

PRECISION/EARDWARE  •  SINGLE  AND  OOU3LE/H32 

.  SINGLE/H36,H48,H60 

REQO.  IMSL  ROUTINES  -  EBALAF,EBBCRF,EHBCKF,EHESSF,EQRE3F,UERTST, 

UGETIO 

NOTATION  >  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 


EIGRF  comoutes  eigenvalues  and  (optionally)  eigenvectors  of  a  real 
natrix.  It  can  also  compute  a  performance  index. 


EIGRF  calls  IMSL  routine  EBALAF  to  balance  th.e  matrix.  Then,  EHESSF 
and  SQRH3F  are  called  to  compute  eigenvalues  and  (optionally)  eigen¬ 
vectors.  When  eigenvectors  are  computed,  EHBCXF  and  EBBCRP'are  called 
to  backtramsform  the  eigenvectors. 


The  performance  index  is  defined  as  follows 


P 


max 


\  \z^\\^  10(N)  (EPS) 


where  the  max  ts  taken  over  the  j  eigenvalues  w.  and  associated  eigen- 

4  J 

.vectors  z-* .  EPS  specifies  the  relative  precision  of  floating  point 
arithmetic .  When  P  is  less  than  1,  the  perfo-rmance  of  the  routines  is 
considered  to  be  excellent  in  the  sense  that  the  residuals  Az-wz  are  as 
small  as  can  be  expected.  When  P  is  between  1  and  100  the  performance 
is  good.  When  P  is  greater  than  100  the  performance  is  considered  poor. 

The  performance  index  was  first  developed  and  used'  by  the  EISPACX  pro¬ 
ject  at  Argonne  National  Laboratory. 
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S««  references: 


1.  Wilkinson,  J.  H.,  The  Algebraic  Eigenvalue  Problem,  Clarendon 
Press,  Oxford,  1963. 


2.  Smith,  B.  T. ,  Boyle,  J.  M.,  Garbow,  B.  3.,  Ikebe,  Y. ,  Klema,  V,  C. 
and  Moler,  C;,  B.,  Matrix  Eigensystem  Routines,  Springer-Verlag, 
1974. 


Programming  Notes 

1.  A  is  preserved  when  IJ0B«2  or  3  2uid  N»IA.  When  N<IA,  rows 
N^*!,  N-1'2,  ...,  lA  of  A  are  destroyed.  In  all  other  cases  A 
is  destroyed. 

e 

2.  The  eigenveilues  are  unordered  except  that  complex  conjugate  pairs 
of  eigenvalues  appear  consecutively  with  the  eigenvalue  having  the 
positive  imaginary  part  first. 

3.  The  eigenvectors  are  not  normalized. 

4.  When  IJ0B«3  (i.e.,  to  compute  a  performance  index  only)  the  eigen¬ 
values,  W,  and  eigenvectors,  Z,  are  assxmed  to  be  input. 

Example  1 


Zn  this  example,  EIGHT  is  called  to  compute  eigenvalues,  eigenvectors 
and  a  perform^Lnce  index  by  setting  input  IJ0B«2.  After  the  call  to 
EZGRF,  the  eigenvectors  are  normalized.  For  machines  which  require 
equivaleneing,  see  example  2. 

INTEGER  N,IA,IJOB,IZ,IER 
HEAL  A(4,4}  ,WK(24) 

COMPLEX  W(4}  ,Z(4,4}  ,ZN 


Input: 


lA 

IZ 

N 

UOB 


A 


4 

4 

4 

2 


4.0 

0.0 

5.0 

3.0 


-5.0  0.0  3.0 

4.0  -3.0  -5.0 

-3.0  4.0  B.O 

0.0  5.0  4.0 


CALL  EIGHT  (A,N,  lA,  IJOB,  W,  Z,  IZ,WE,  lER) 

NORMALIZE  EIGENVECTORS 


DO  5  J«1,N 
ZN  «  Z(I,J) 

DO  5  I»1,N 

Z(Z,J)  »  Z(I,J)/ZN 
)  CONTINUE 
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Output 
lER  ■ 


0 


W 

Z  * 


( 12. 0/1. 0't‘5 101, 1.0~5. 01,2.0)  ( eigenvalues ) 


1.0 

1.0 

0.0 

1.0 

0.0 

1.0 

-1.0 

0.0 

+1 

-1.0 

0.0 

+1 

1.0 

1.0 

1.0 

9 

0.0 

-1.0 

9 

0.0 

1.0 

9 

-1.0 

1.0 

-1.0 

0.0 

-1.0 

0.0 

1.0 

(eigenvectors) 


WK(1)  <  1.0  (performance  Index) 

Example  2 

For  machines  which  require  an  equivalence  statement  In  situations  where 
an  array  Is  of  one  type  In  the  calling  program  but  of  another  type  In 
the  subroutine,  EIGRF  should  be  called  as  follows. 


INTEGER 

REAL 

COMPLEX 

EQUIVALENCE 


N,IA,IJOB,IZ,IER 

A(4,4) ,WK(24) ,RW(8) ,RZ(32) 

W{4) ,2(4,4) ,2N 

(W(l)  ,RW(1))  ,  (2(1,1)  ,RZ(1)) 


Input 

• 

• 

lA 

« 

4 

IZ 

m 

4 

N 

m 

4 

UOB 

m 

2 

-5.0  0.0 

• 

4.0 

% 

0.0 

4.0  -3.0 

A 

5.0 

-3.0  4.0 

3.0 

0.0  5.0 

CALL 

EIGRF  (A,N, 

IA,IJOB,RW,RZ 

C 

00  5 

J-1 

,N 

ZN 

a 

Z(1,J) 

DO 

5 

I« 

>1,N 

Z(I,J)  >  2(I,J)/2N 
5  CONTINUE 


3.0 

-5.0 

0.0 

4.0 

I2,WK,ISR) 

NORMALIZE  EIGENVECTORS 


Output: 
ZZR  •  0 


W  •  (12.0,1.0-^5.01,1.0-5.01,2.0)  (eigenvalues) 


1.0 

1.0 

0.0 

1.0 

0.0 

1.0 

-1.0 

0.0 

+1 

-1.0 

0.0 

+1 

1.0 

1.0 

1.0 

9 

0.0 

-1.0 

9 

0.0 

1.0 

9 

-1.0 

»-• 

. 

o 

1 _ 

-1.0 

0.0 

-1.0 

0.0 

1.0 

NKd)  <  1.0  (performance  Index) 


(eigenvectors) 
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IMSL  ROUTINE  NAME 


LEQT2F 


PURPOSE  -  LINEAR  EQUATION  SOLUTION  -  FULL  STORAGE 

MODE  -  HIGH  ACCURACY  SOLUTION 

USAGE  -  CALL  LEQT2F  (A,M,N,IA,B,IDGT,WKAREA,IER) 

ARGUME?!TS  A  -  INPUT  MATRIX  OP  DIMENSION  N  BY  N  CONTAINING 

THE,  COEFFICIENT  MATRIX  OF  THE  EQUATION 
AX  a  B. 

M  '  -  NUMBER  OF  RIGHT-HAND  SIDES.  (INPUT) 

N  -  ORDER  OF  A  AND  NUMBER  OF  ROWS  IN  B.  (INPUT) 

lA  -  ROW  DIMENSION  OF  A  AND  B  EXACTLY  AS  SPECIFIED 

IN  THE  DIMENSION  STATEMENT  IN  THE  CALLING 
PROGRAM.  (INPUT) 

B  .  -  INPUT  MATRIX  OF  DIMENSION  N  BY  M  CONTAINING  ■ 
THE  RIGHT-HAND  SIDES  OF  THE  EQUATION  AX  =  B 
ON  OUTPUT,  THE  N  BY  M  MATRIX  OP  SOLUTIONS 
REPLACES  B. 

IDGT  -  INPUT  OPTION. 

IF  IDGT  IS  GREATER  THAN  0,  THE  ELEMENTS  OF 
A  AND  B  ARE  ASSUMED  TO  BE  CORRECT  TO  IDGT 
DECIMAL  DIGITS  AND  THE  ROUTINE  PERFORMS 
AN  ACCURACY  TEST. 

IF  IDGT  EQUALS  0,  THE  ACCURACY  TEST  IS 
BYPASSED. 

ON  OUTPUT,  IDGT  CONTAINS  THE  APPROXIMATE 
NUMBER  OF  DIGITS  IN  THE  ANSWER  WHICH 
WERE  UNCHANGED  AFTER  IMPROVEMENT. 

WKAREA  -  WORK  AREA  OF  DIMENSION  GREATER  THAN  OR  EQUAL 
TO  N**2+3N. 

lER  -  ERROR  PARAMETER.  (OUTPUT) 

WARNING  ERROR 

,  lER  «  34  INDICATES  THAT  THE  ACCURACY  TEST 
FAILED.  THE  COMPUTED  SOLUTION  MAY  BE  IN 
ERROR  BY  MORE  THAN  CAN  BE  ACCOUNTED  FOR 
BY  THE  UNCERTAINTY  OF  THE  DATA.  THIS 
WARNING  CAN  BE  PRODUCED  ONLY  IF  IDGT  IS 
GREATER  THAN  0  ON  INPUT.  (SEE  THE 
CHAPTER  L  PRELUDE  FOR  FURTHER  DISCUSSION. 
TERMINAL  ERROR 

lER  »  129  INDICATES  THAT  THE  MATRIX  IS 
ALGORITHMICALLY.  SINGULAR.  (SEE  THE 
CHAPTER  L  PRELUDE) . 

lER  «  131  INDICATES  THAT*  THE  MATRIX  IS  TOO 
ILL-CONDITIONED  FOR  ITERATIVE  IMPROVEMENT 
TO  BE  EFFECTIVE. 

PRECIS lON/HARDWARE  -  SINGLE  AND  DOUBLE/H32 

-  SINGLE/H36,H48,H60 

REQD.  IMSL  ROUTINES  -  SINGLE/LUDATN,LUELMN,LUREFN,UERTST,UGETIO  ‘ 

-  DOUBLE/LUDATN , LUELMN , LUREFN , UERTST , UGETIO , 

VXADD , VXMUL , VXSTO 

NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

June,  1982 
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LEQT2F-1 


Algorithm 

LEQT2F  solves  the  set  of  linear  equations  AX^B  for  X,  where  A  is  the 
N  by  N  matrix  and  is  in  full  storage  mode.  B  is  N  by  M.  The  different, 
between  this  routine  and  routine  LEQTIF  is  that  LEQT2F  invokes  iterative 
improvement  if  necessary,  in  order  to  improve  the  accuracy  of  the 
solution  X. 

The  routine  performs  Gaussian  elimination  (Grout  algorithm)  wifth  equi¬ 
libration,  partial  pivoting,  and  iterative  improvement  as  required. 

See  reference: 

Forsythe,  George  and  Mcler,  Cleve  B.,  Computer  Solution  of  Linear 
Mgebraic  Systems,  Englewood  Cliffs,  N.  J. ,  Prentice-Hall ,  Inc.,  1967, 
Chapters  9,  13,  24. 

Programming  Notes 

1.  Iterative  improvement  is  costly  in  both  computer  time  and  storage. 
When  high  accuracy  is  not  needed,  subroutine  LEQTIF  may  be  used  to 
advantage. 

2.  When  lA  is  greater  than  N,  elements  of  A  in  rows  N-(>1  to  lA  are 
used  as  workspace  and  are  destroyed.  However,  the  first  N  rows 
of  A  are  restored  to  their  original  content  on  exit  from  LEQT2F. 

Accuracy 

If  I06T  is  greater  than  zero,  elements  of  A  are  assumed  to  be  correct 
to  IDGT  decimal  digits.  The  solution  X  will  be  the  exact  solution, 
without  any  roundoff  error,  to  a  matrix  X  whose  elemehts  agree  with 
the  elements  of  A  in  the  first  IDGT  decimal  digits.  The  program  first 
attempts  such  a  solution  without  iterative  improvement.  Then  iterative 
improvement  is  performed  if  necessary.  If  this  also  fails,  solution  is 
not  possible  ^md  the  program  exits.  Upon  exit,  the  first  columns  of  B 
will  have  been  replaced  by  the  best  solution  that  the  computer  can 
generate  and  IDGT  is  set  to  the  approximate  number  of  digits  in  the . 
emswer  which  were  unchanged  by  the  improvement  (see  IMSL  routine  LUREFF) 
The  other  columns  of  B  are  left  unchanged  in  this  case  and  lER  is  set  to 
131.  If  input  IDGT  equals  zero,  iterative  improvement  is  automatically 
perfozmed. 

Examp le 

This  example  inputs  the  3  by  3  matrix  A  and  the  3  by  4  matrix  B  solving 
for  the  3  by  4  matrix  X  of  AX«B.  X  overwrites  B  on  output. 

Input: 

REAL  A(4,4),B(4,4),WKAR£A(18) 

INTEGER  M,N,IA,IOGT,IER 
N  «  3 

M  «  4 

lA  -  4 

IDGT  -  3 
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A 


33.000 

16.0 

72.0  X 

- 

-24.000 

-10.0 

57.0  X 

-  8.000 

-  4.0 

17.0  X 

1 

X 

X 

X  X 

1.0 

0.0  0.0 

-359.0 

0.0 

1.0  0.0 

281.0 

0.0 

0.0  1.0 

85.0 

X 

X  X 

X 

CALL  L£QT2F(A,M,N, IA«B, IDGT,WKAR£A,IER) 


Output: 

IDGT  = 

3 

e 

lER  - 

3 

-9.66666 

-2. 66667 

-32. 

1. 

B  - 

8.0 

2.  5 

25.5 

-2. 

2.66667 

.666667 

9. 

-5. 

X 

X 

X 

X 

Note:  X  indicates  elements  not  used  by  LEQT2F. 


P 
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IMSL  ROUTINE  NAME 


EIGRS 


PURPOSE 

USAGE 

ARGUMENTS 


-  EIGENVALUES  AND  (OPTIONALLY)  EIGENVECTORS  OF 

A  REAL  SYMMETRIC  MATRIX 

-  CALL  EIGRS  (A,N, JOBN,D,Z,IZ,WK,IER) 

A  -  INPUT  REAL  SYMMETRIC  MATRIX  OF  ORDER  N, 

WHOSE  EIGENVALUES  AND  EIGENVECTORS 
ARE  TO  BE  COMPUTED.  INPUT  A  IS 
DESTROYED  IF  IJOB  IS  EQUAL  TO  0  OR  1 . 

N  -  INPUT  ORDER  OF  THE  MATRIX  A. 

JOBN  -  INPUT  OPTION  PARAMETER.  IF  JOBN.GE.IO 
A  IS  ASSUMED  TO  BE  IN  FULL  STORAGE  MODE 
(IN  THIS  CASE,  A  MUST  BE  DIMENSIONED  EXACTLY 
,  N  BY  N  IN  THE  CALLING  PROGRAM) . 

IF  JOBN. LT. 10  THEN  A  IS  ASSUMED  TO  BE  IN 
SYMMETRIC  STORAGE  MODE.  DEFINE 
IJOB=MOD(JOBN,10) .  THEN  WHEN 

IJOB  =  0,  COMPUTE  EIGENVALUES  ONLY 
IJOB  =  1,  COMPUTE  EIGENVALUES  AND  EIGEN¬ 
VECTORS  . 

IJOB  =  2,  COMPUTE  EIGENVALUES,  EIGENVECTORS 
AND  PERFORMANCE  INDEX. 

IJOB  =  3,  COMPUTE  PERFORMANCE  INDEX  ONLY. 

IF  THE  PERFORMANCE  INDEX  IS  COMPUTED,  IT  IS 
RETURNED  IN  WK(1).  THE  ROUTINES  HAVE 
PERFORMED  (WELL,  SATISFACTORILY,  POORLY)  IF 
WK(1)  IS  (LESS  THAN  1,  BETWEEN  1  AND  100, 
GREATER  THAN  100) . 

D  -  OUTPUT  VECTOR  OF  LENGTH  N, 

CONTAINING  THE  EIGENVALUES  OF  A  IN  ASCENDING 
ORDER. 

Z  -  OUTPUT  N  BY  N  MATRIX  CONTAINING 

THE  EIGENVECTORS  OF  A. 

THE  EIGENVECTOR  IN  COLUMN  J  OF  Z  CORRES¬ 
PONDS  TO  THE  EIGENVALUE  D(J). 

IF  IJOB  =0,  Z  IS  NOT  USED. 

IZ  -  INPUT  ROW  DIMENSION  OF  MATRIX  Z  EXACTLY  AS 

SPECIFIED  IN  THE  DIMENSION  STATEMENT  IN  THE 
CALLING  PROGRAM. 

WK  -  WORK  AREA,  THE  LENGTH  OF  WK  DEPENDS 

ON  THE  VALUE  OF  IJOB,  WHEN 

IJOB  *  0,  THE  LENGTH  OF  WK  IS  AT  LEAST  N. 

IJOB  =  1,  THE  LENGTH  OF  WK  IS  AT  LEAST  N. 

IJOB  =  2,  THE  LENGTH  OF  WK  IS  AT  LEAST 

N(N+l)/2+N. 

IJOB  =  3,  THE  LENGTH  OF  WK  IS  AT  LEAST  1. 
lER  -  ERROR  PARAMETER  (OUTPUT) 

TERMINAL  ERROR 

lER  «  12 8+ J,  INDICATES  THAT  EQRT2S  FAILED 
TO  CONVERGE  ON  EIGENVALUE  J.  EIGENVALUES 
AND  EIGENVECTORS  1,...,J-1  HAVE  BEEN 
COMPUTED  CORRECTLY,  BUT  THE  EIGENVALUES 
ARE  UNORDERED.  THE  PERFORMANCE  INDEX 
IS  SET  TO  1000.0 


November,  1984 
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L-  PRECISION/HARDWARE 


REQD.  IMSL  ROUTINES 
NOTATION 


WARNING  ERROR  (WITH  FIX) 

IN  THE  FOLLOWING,  IJOB  =  MOD (JOHN, 10) . 
lER  =  66,  INDICATES  IJOB  IS  LESS  THAN  0  OR 
IJOB  IS  GREATER  THAN  3.  IJOB  SET  TO  1. 
lER  =  67,  INDICATES  IJOB  IS  NOT  EQUAL  TO 
ZERO,  AND  IZ  IS  LESS  THAN  THE  ORDER  OF 
MATRIX  A.  IJOB  IS  SET  TO  ZERO. 

SINGLE  AND  DOUBLE/H32 
S INGLE/H  36,H48,H60 

EHOBKS , EHOUSS , EQRT2S , UERTST , UGETIO 

INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUJZTION  OR  THROUGH  IMSL  ROUTINE  UHELP 


Algorithm 

EIGRS  calls  IMSL  routine  EHOUSS  and  EQRT2S  to  compute  eigenvalues  and 
(optionally)  eigenvectors.  When  eigenvectors  are  computed,  EHOBKS  is 
called  to  backtransform  the  eigenvectors. 

The  performance  index  is  defined  as  follows: 


P  = 


max 

l£j^n 


|Az^-d^2^ I  1^ 


10(N) (EPS) 


te  the  max  is  taken  over  the  n  eigenvalues  d^  and  associated  eigen¬ 
vectors  z^.  EPS  specifies  the  relative  precision  of  floating  point 
arithmetic.  When  P  is  less  than  1,  the  performance  of  the  routines  is 
/  considered  to  be  excellent  in  the  sense  that  the  residuals  Az-dz  are  as 
'  small  as  can  be  expected.  When  P  is  between  1  and  100  the  performance 
I  is  good.  When  P  is  greater  than  100  the  performance  is  considered  poor, 

The  performance  index  was  first  developed  and  used  by  the  EISPACK  pro- 
ject  at  Argonne  National  Laboratory. 

■  See  references: 


1. 


12. 


Wilkinson,  J.  H.,  The  Algebraic  Eigenvalue  Problem,  Clarendon  Press 
Oxford,  1965. 

Smith,  B.  T.,  Boyle,  J.  M. ,  Garbow,  B.  S.,  Ikebe,  Y.,  Klema,  V.  C., 
and  Moler,  C.  B.,  Matrix  Eigensystem  Routines,  Springer-Verlag, 
1974. 


Programniing  Notes 

1.  A  is  preserved  when  IJOB=2  or  3  (IJOB=MOD(JOBN,10) ) .  In  all 
other  cases  A  is  destroyed. 

2^.  The  computed  eigenvectors  are  normalized  to  each  have  Euclidean 
^ ■  length  1 . 
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IMTEGEH  N, JOBN, IZ, lER 

REAL  A(10)  ,0(4)  ,Z(4,  4)  ,WK(14} 

N  «  4 

IZ  «  4 


A 


5.0 

4.0  5.0 
1.0  1.0 
1.0  1.0 


4.0 

2.0 


In  symmetric  storage  mode  A(l)>5.0, 
A(2)-4.0,A(3)-5.0, ...,A(10)-4.0 


JOBN  «  2 

CALL  EIGRS  (A,  2T,  JOBN,  0,  Z,  IZ*,  WK,  lER} 


Output: 


lER  «  0 

0  *  (1.0,  2.0,  5.0,  10.0)  (eigenvalues) 


Z 


0.70710 

0.70710 

0.00000 

0.00000 


0.00000 

0.00000 

0.70710 

-0.70710 


-0.31622 

-0.31622 

0.63245 

0.63243 


0.63245 

0.63245 

0.31622 

0.31622 


(eigenvectors) 


WR(1)  <1  (performance  index) 

Note:  Z  is  uniq^io  to  within  a  sign  change  for  each  column. 
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tlng  their  efforts  with  the  method.;  The  eigenvalue  method  Is  used 
to  solve  the  parabolic  heat  equation  In  multidimenslons  with  .  „  i. 

transient  temperatures.  Extensions,  Into  three  dimensions  are  made 
to  determine  the  method's  feasibility  In  handling  large  geometry 
problems  requiring  great  numbers  of  Internal  mesh  points. 

The  eigenvalue  method  proves  to  be  slightly  better  In  accuracy 
than  the  comparison  routines  because  of  an  exact  treatment,  as 
opposed  to  a  numerical  approximation,  of  the  time  derivative  In  the 
heat  equation.  It  Is  an  unconditionally  stabel  method.  It  has  the 
potential  of  being  a  very  powerful  routine  In  solving  long  transient 
type  problems.  The  method  Is  not  well  suited  to  finely  meshed  grid 
arrays  or  large  regions  because  of  the  time  and  memory  requirements  , 
necessary  for  calculating  large  sets  of  eigenvalues  and  eigenvectors. 
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